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Abstract 



We present an introduction to the study of chaos in discrete and continuous dynam- 

i 1 ical systems using the CAS Maxima. These notes are intended to cover the standard 

Q topics and techniques: discrete and continuous logistic equation to model growth popu- 

lation, staircase plots, bifurcation diagrams and chaos transition, nonlinear continuous 
dynamics (Lorentz system and Duffing oscillator), Lyapunov exponents, Poincare sec- 
tions, fractal dimension and strange attractors. The distinctive feature here is the use 
of free software with just one ingredient: the CAS Maxima. It is cross-platform and 
have extensive on-line documentation. 
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1 Introduction 



The use of a computer is a natural way to explore those notions related to dynamical 
systems; here, we present the basics of chaotic dynamics using Maxima. This is a free, 
cross-platform, general purpose computer algebra system |Max capable also of doing 
numerical computations, that we have used in short introductory summer courses (about 
10 hours). The following notes are the outcomes of such courses, and given the target 
audience, we have made every effort to keep them at an elementary level, relegating the 
most technical points to the footnotes. 

Of course, it is impossible to include all topics one would wish, so we have excluded 
those theoretical results which are hard to implement using the computer, because they 
involve irrational numbers. Explicitly, we only make some brief remarks about results such 
as the density of periodic orbits or sensitivity to initial conditions and leave aside topological 
transitivity. 

In writing this paper, we have used Maxima version 5.24.0 and its graphical interface 
wxMaxima version 11.04.0. The version of Gnuplot was 4.4 patchlevel 3, everything running 
on a generic desktop PC with Slackware Linux 13.1. 



2 The logistic equation: continuous case 

The simplest model of population growth is given by: 

dP(t) 



dt 



k ■ P(t) (1) 



where P(t) is the population at time t and k is a constant, positive for an increasing 
population and negative for a decreasing one. However, in the long run this is not a good 
model: it disregards limiting factors such as propagations of diseases or lacks of food supply. 
A simple modification in |l]) which takes into account these factors can be obtained replacing 
k with a parameter K = K(P), such that it decreases when P increases. A possible model 
for this, is K = a — b ■ P, a S> b > 0, so ([I]) takes the form 

d ^ = (a-b-P(t))-P(t). (2) 

If P{t) ~ 0, equation Q becomes equation Q, with k ~ a, but for increasing values of 

therefore dP ^ 

Making the change of variables 



P(t), b ■ Pit) approach a, and therefore d ^ ' ~ 0, slowing down the population's growth. 



ax dP(t) adx 



equation Q becomes (with a = k) 

n t 

- = k-x-(l-x), (4) 

which is called the logistic equation (it was introduced by P. F. Verhulst in 1838, a discrete 
version was popularized by R. May in the seventies, see [May 76| ). 
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We will use Maxima to solve equation Q (note the use of the apostrophe to define an 
equation. The command for solving first or second-order ODEs is ode2): 

C/„il) Miff (x,t)=k*x*(l-x) ; 

(%ol) — x = k (1 — x) x 
(%i2) ode2C/.,x,t); 

(%o2) iog(x) -;° g(x - i) =t + %c 

k 

Now, we set an initial condition for this first-order (hence the 1 in icl) equation x(t = 0) = 
x : 

(%i3) icl(7„,t=0,x=x[0] ) ,logcontract; 
log(^) kt + logt^- 



(%o3) 



XQ — 1 



k k 
and solving the last equation we get 
(%i4) solve (%, x) ; 

(%° 4 ) [* = ^^^J 
xo e K 1 — xo + 1 

Thus, we have the explicit solution (note that the outcome (%o4) is a list -a collection of 
elements enclosed between brackets- and we want to define our solution as the right-hand 
side of its first element): 

(%i5) define(x(t) ,rhs(f irst (%o4) ) )$ 

Let us note that normalization ^ implies that for all initial conditions < xq, the 
solutions x(t) — > 1 when t —> oo. In fact, even if we had not known the solutions explicitly, a 
simple continuity argument would show that all solutions bounded by the steady states x = 
and x = 1, tend monotonically to x = Also, in the unbounded region x > 1 all solutions 
tend monotonically to the steady state x = 1. For x < there is no solution (populations 
are non negative). In general, in a one dimensional continuous system, solutions either 
converge monotonically to a steady state or diverge. 

The parameter k affects the slope of the solutions, as we can see by graphing them 
together for different values of k (in this example, varying from k = 0.25 to k = 1.75 with 
step 0.25). 



x For a system x — f(x), in the region between two consecutive steady states x(t) — x* with f(x*) = 0, 
x* < x%, the mapping / (and therefore the derivative x) has constant sign, so a solution x(t) with initial 
condition xo such that x\ < xo < x\, is strictly monotonic. This solution is bounded, so linit-^oo x(t) = x% if 
x(t) is increasing, and limt-joo x(t) — x* if x(t) is decreasing. Finally, because of the uniqueness theorem for 
the Cauchy problem, x(i) can not reach the values x\,X2- A similar reasoning can be used for unbounded 
regions. 
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(%i6) wxplot2d(makelist(subst( [k=d*0.25,x[0]=0. 1] ,x(t)) ,d,l,7) , 
[t,0,15] , cons (legend, makelist (k=d*0 . 25, d, 1,7)) , 
[gnuplot_preamble , "set key right bottom"]); 




2 4 6 8 10 12 14 



(%t6) * (%o6) 

3 The logistic equation: discrete case 

The discrete version of equation Q is the difference equation (with r > 0) 

x n+ i = r ■ x n ■ (1 - x n ). (5) 

In general, equation ^ can't be solved exactljj^J Let us define the so called logistic map 
(in fact, a family of quadratic maps), F(r,x) = r ■ x • (1 — x). From §5§ we see that, for 
a given r, the fixed points of F(r,x) are the only constant solutions x n = c of d5J). These 
steady states can be easily calculated: 

(%i7) solve (r*c* (l-c)=c , c) ; 

(%o7) [ c = 1 —±,c = Q] 

From our experience with the continuous case, we are tempted to say that any solution 
of the logistic map is either a sequence bounded by the constant solutions x n = and 
x n = (r — l)/r (asymptotically approaching one of these), or a divergent one. But this is 
not true: The behavior of discrete systems is quite different from that of their continuous 
counterparts. To be sure, we make a graphical analysis of the solutions. First, we define 
the evolution operator of the system: 

(°/.i8) F (r , x) : =r*x* ( 1-x) $ 

Next, we use the evolution command (included in the dynamics package) to calculate, 
for a fixed r, the n + 1 points Xj+i = F(r,Xi) from i = to i = n, where the initial value 
xq is given (here we use a pseudo-random initial condition between and 1). 

2 This sentence deserves a clarification. Actually, "closed" solutions do exist, although defined throuh 
implicit functional formulas which are 'hardly useful for computational purposes', see IBru 98) and ref- 
erences therein. However, there are well-known explicit solutions for particular values of k, see http: 
//en. wikipedia. org/wiki/Logistic_map#Solution_in_some_cases 
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This set of points is a segment of the orbit of xq. We start with r = 0.25 making 15 
iterations: 

(%i9) load (dynamics) $ 

(%ilO) set_random_state(make_random_state (true))$ 
('/.ill) x[0] : random (1.0) ; 
(%oll) .8886589561406515 

C/.il2) evolution(F(0.25,x) ,x[0] ,15, [y,0,l]) ; 



%tl2) 




10 12 14 



(%ol2) 



In this case, the population dies. Indeed, this unlucky fate is independent of the initial 
condition xq , whenever < r < 1 and xq G ] r - = ^ , £ [ (let us note that this last condition 
is always fulfilled if xq g]0,1[). For < r < 1, this parameter controls the rate at which 
population dies, as the reader can check by using other values for r, such as r = 0.123456789. 

A few experiments will lead us to the following rules for the behavior of the solutions: 

• If < r < 1 and xq G ] £ [, the population eventually dies, independently of the 
initial condition. 

• If 1 < r < 2 and xo G]0, 1[, the population quickly approaches the value (r — 1)/V, 
independently of the initial condition. 

• If 2 < r < 3 and xq e]0, 1[, the population tends again to the value (r — l)/r, but in 
an oscillating way (maybe after a short transient). The rate of convergence is linear, 
except for r = 3, where the rate of convergence is quite slow, in fact it is sub-linear. 

• If3<r<l + \/6 (with 1 + \/6 ^ 3.45), the population oscillates between two values, 
almost independently of the initial condition. These two values, which depend on r, 
are said to have primary period two. 

• If 3.45 < r < 3.54 (approximately), for almost all initial condition the population 
oscillates between four periodic points. 

• If r increases over 3.54, for almost all initial condition the population oscillates between 
8 periodic points, then 16, 32, etc. The size of the intervals formed by the values of the 
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parameter producing oscillations of a given length, becomes small, and the quotient 
of the size for two consecutive period-doubling intervals approaches the so called 
Feigenbaum constant F = 4.669... This behavior is called a period- doubling cascade. 

Let us note that when the parameter r increases, the dynamics differs from that observed 
in the continuous case. 

Let us show some examples. For our initial pseudo-random xq and r = 1.3: 

(°/.il3) evolution(F(1.3,x) ,x[0] ,15, [y,0,l]) ; 



1 

G.8 
0.6 
0.4 
0.2 





2 4 6 



10 12 14 



(%tl3) 



%ol3) 



We note that (1.3 — 1)/1.3 ~ 0.23. For r = 3.5 and the initial condition xq = 0.3, we 
observe the transient and four periodic points: 

(°/„il4) evolution(F(3.5,x) ,0.3,25, [y,0,l]) ; 



1 

0.8 
0.6 
0.4 
0.2 




• • • ♦ 
• • • • • 



(%tl4) 



5 10 15 20 25 

n 



(%ol4) 



4 Staircase diagrams 

There is another graphical method to study the phenomenon of the emergence of "attract- 
ing" periodic points that we have seen in the last example, based on the particular form of 
the evolution equation: x n+ \ = f(x n ). For the logistic map (f(x) = r ■ x ■ (1 — x)), if the 
coordinates of the points in the plane represent two consecutive values in the orbit of xq, 
the point (x,y) = (x n ,x n +i) can be obtained graphically as the intersection of the vertical 
line x = x n and the graph of the function y = f(x). Once x n +i is known, in order to get 
x n +2 we just have to let play the role of x n in the previous step. So, we take the 
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intersection of the horizontal line y = with the diagonal y = x. Now x n+ 2 = f(x n +\), 
and by iterating the process a predetermined number of times, we get a segment of the orbit 
of Xq. Maxima implements this construction with the staircase command: 

(°/ il5) staircase (F (3. l,x) ,x[0] ,250, [x,0, 1] , [y,0, 1] ) ; 




0.2 0.4 0.6 0.8 1 



(%t!5) x(n) (%ol5) 



Here we see the orbit oscillating between the periodic points (x ~ 0.56) and (x ~ 0.76). 
We can also reproduce the behavior of (%ol4), where there is a period- four orbit (we make 
500 iterations to get a distinguishable path). Let us note that each "rectangle" intersecting 
the diagonal defines two periodic points: 

(%il6) staircase(F(3.5,x) ,x[0] ,500, [x,0,l] , [y,0,l]) ; 



%tl6) 




(%ol6) 



5 Bifurcations and chaos 

We have seen that the behavior of the orbits for the logistic map depends on the parameter 
r. In fact, from r = 3 onwards, the number of periodic points increases in a kind of "period- 
doubling cascade" (see (%ol7) below). It is possible to study the change in the structure 
of the periodic orbits using the so called bifurcation diagrams. In these, we represent the 
values of the parameter r on the horizontal axis. For each of these values we mark the values 
of the corresponding "attracting" periodic points (the use of attracting points is a technical 
issue, that we will not address here). Maxima implements these plots in its built-in orbits 
command. In the example below, we restrict the plot to the iterations between n = 150 
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and n = 200, for a given initial condition xq. The range of the values of r is [2.5,4]: 



(%il7) orbits (F(r,x) , x[0], 150, 200, [r, 2.5, 4], [style, dots]); 



%tl7) 
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(%ol7) 



The diagram so obtained displays a special property: it is self-similar. This means that 
after a change of scale, the resulting picture has the same structure that the original one. 
For example, if we magnify (%ol7) centering the zoom on the second bifurcation of the 
lower branch, we get: 

(7.il8) orbits(F(r,x) ,x[0] ,150,200, [r,3.5,3.6] , [x, 0.3, 0.4] , 
[style, dots] ) ; 



0.4 
0.38 
36 
0.34 
0.32 



03 I ■ ■ ■ ■ 1 

3 5 3 52 3 54 3 56 3 58 3 6 

(%ol8) 



%tl8 



A measure of the self-similarity of a set is its fractal dimension which, unlike the dimen- 



sion of a vector space, can be a non integer number. In Section 10 we show an example of 
a set having non integer fractal dimension. 

The period-doubling cascade and the appearance of sets with fractal dimension, are 
usually signs of chaotic behavior. A dynamical system with whose evolution is described 
by a continuous function / : I C K — > I is said to be c/iaoiicj^] if: 

(a) It has sensitivity to initial conditions: a small variation in the initial condition xq of an 
orbit can produce, in the long run, another orbit far away from the original one. 



3 For simplicity we consider only one-dimensional systems, but the case of multidimensional systems is 
analogous. For further details on the chaos definition see Ba n 92] . [Dev 89j . |VB 94] 
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(b) Periodic points are "dense" in the phase space I of the system: in general, a set S is 
dense in the set I C K if for any point p G I and any e > 0, the interval ]j> — e,p + e[ 
intersects S. 

(c) It has the mixing property (we also say that the system is topologically transitive): for 
any two intervals J,K C I there exist points of J whose orbits eventually enter K, that 
is, there exists an n > 1 such that 



n) 



with / n = / o • • • o /. This property is satisfied if and only if the system has an orbit 
{fn(x ) : n £ N U {0}} that is dense in I. 

We illustrate each one of these conditions using the logistic map in the examples below. 



(a) Sensitivity to initial conditions. For r = 4, we consider the initial conditions xq 
= | and xq + e, comparing the respective orbits after 50 iterations. 

(%il9) evolution(F(4,x) ,3/4,50, [style, [lines, 2]] ) ; 



758 
0.756 
0.754 
0,752 
0.75 
0.748 
0.746 
0.744 
742 



50 



10 20 30 40 

(%tl9) " (%ol9) 

Choosing e small enough (pseudo-random) we see how drastically these orbits diverge 
(we will see how to measure this divergence in section l7j): 

(%i20) eps: random (0. 000000000001) ; 
(%o20) 9.87415729237455 1(T 13 

(%i21) evolution(F(4,x) ,3/4+eps,50, [style, [lines, 2]] ) ; 



4 The equivalence is as follows: if the orbit of xq is dense in J, there exists an infinity of points x n = f n {xo) 
in any pair of intervals J, K, so / is transitive. Conversely, if / is transitive, for a given interval J let {Aj}j^fs 
be the basis of open intervals with rational end points that are contained in J; from the transitivity of / 
the family of open sets Bj = U P £E,f- P (A n ) is dense in J (since, if U C J, there exists a number p such that 
f p (U)(~}A n / 0, that is, U n f p (A n ) = UC]B n ^ 0)- From the Baire theorem, (IjesqBj is dense in J too. But 
this intersection consists of points with dense orbits, since if x G Hj^Bj, for each n S N there is a p n G N 
such that f Pn £ A„. The sequence {/ Pn }„a has points in each A n and therefore is dense. 
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(b) The set of periodic points is dense in the interval / = [0,1]. This is clearly seen in 
(%ol7), where the right side of the unit square is "filled" with periodic points. 

(c) There exists a dense orbit. Again, we can exemplify this behavior with k = 4 and the 
initial condition xq = 0.2 after 50, 000 iterations. Its orbit fills the phase space: 



(%i22) evolution(F(4,x) ,0.2,50000, [style, dots] ) ; 



0.8 

0.7 




D 10000 20000 30000 40000 50000 



(%o22) 



(%t22) 



6 The Lorenz attractor 

As mentioned, the monotonicity of the solutions restricts the presence of chaos in one- 
dimensional continuous systems. However, for higher order equations (at least of second 
order with a time dependent non-homogeneous term, as in the Duffing equation shown 
later) or in systems of first order equations (at least three equations, even in the case of 
an autonomous system, but necessarily non lineaij^]), it is possible to have chaotic behavior. 
As an example, we consider the Lorenz system, the first system in which the properties of 



5 We have already seen that one dimensional continuous systems can not display chaotic behavior (nor 
even periodic behavior). For autonomous systems in the plane (or second order differential equations with 
time-independent coefficients), the uniqueness of solutions together with the Poincare-Bendixson theorem 
rule out the existence of chaos. The problem of setting the "minimal conditions" to ensure the existence of 
chaos is still an open question. An interesting report on this topic is |SL 00] , 
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chaotic dynamics were explicitly seen, in 1963 |Lor 63j : 



x = lOy — Wx 
< y = —xz + 28x — y 
k z = xy — 8z/3 

To numerically study the orbits, we use the 4th-order Runge-Kutta method. The following 
commands compute the orbit with initial condition (—8, 8, 27) after 50 time units (with a 
step of 0.01): 

(%i23) numer : f alse$ 

(%i24) latractor: [10*y-10*x, -x*z+28*x-y, x*y-8*z/3] $ 
(7„i25) linitial: [-8, 8, 27]$ 

C/.i26) lsolution:rk (latractor, [x,y,z] , linitial, [t, 0,50,0. 01] )$ 
(%i27) lpoints: map (lambda ( [x] , rest(x)), lsolution)$ 
(°/ i28) load (draw) $ 

The following command needs that the package draw be previously loaded (do load (draw) 
otherwise). Note that the output is actually a separate Gnuplot window, so (by pressing 
the left button of the mouse and dragging around) the reader can rotate the whole picture 
to better appreciate its details. 

(%i29) wxdraw3d(point_type=none , point s_ j oined=true , color=orange , 
xlabel="x(t) " ,ylabel="y(t) " ,zlabel="z(t) " , 
xtics=10 ,ytics=10 ,ztics=10, points (lpoints) ) ; 




(%t29) (%o29) 

The phrase "order within chaos" is frequently found. It means that in some systems, even 
though the orbits can appear unpredictable and very complex, there is some regularity. For 
example, in the Lorenz system, the sensitivity to initial conditions prevents us from making 
accurate predictions, but we can assert that in the long run the orbit will be restricted to a 
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bounded region of the space called an aiiracfor[^] This seems clear in figure (%o29), but we 
can prove it analytically with the aid of Maxima. Let us consider a general Lorentz system, 
where a, p and /3 are parameters: 

x = a(y — x) 

y = x (p - z) - y , (6) 

z = xy — (3z 

(%i30) remfunction(x)$ 

(°/„i31) gradef (x(t) ,%sigma* (y (t) -x(t) ) ) $ 

(°/„i32) gradef (y(t),x(t)*(%rho -z(t) )-y(t) )$ 

(%i33) gradef (z(t) ,x(t) *y(t)-°/ beta *z(t))$ 

We consider the set of those points belonging to an orbit that, for a given instant t, 
lie on the sphere with center (0, 0, a + p) and radius R(a + p), with R a constant to be 
determined: 

(°/„i34) x(t)-2+y(t)-2+(z(t)-(%sigma +%rho) ) ~2=R"2* (%sigma +7,rho)-2; 

(%o34) (z (t) - a - pf + y {tf + x {tf = (a + pf R 2 

Differentiating equation (%o34) (using the chain rule) and substituting ([6]), we get: 
(°/.i35) facsum(diff (°/.,t) ,x(t) ,y(t) ,z(t)) ; 

(%o35) -2f3z(t) 2 + 2(3 (a + p) z (t) - 2y {tf - 2ax (t) 2 = 

This looks like an ellipsoid, and we can prove it is one by "completing the squares": 

(°/ i36) factor (solve ( (sqrt (2*%beta) *z (t) -u) "2-u"2= 

2*7,beta*z(t) "2-2*y o beta*(y,sigma+y o rho)*z(t) ,u)) ; 



(%o36) [u = -j= ] 

So, the points (x(t),y(t), z(t)) should satisfy 



-2ax(t) 2 - 2y(t) 2 - 2/3 ( z(t) - Z±f) + = °' 



which is the equation of an ellipsoid with center (0,0, (cr + p)/2). The semi-axis along 
directions OX, OY and OZ are, respectively: 



6 Attractors are classified as strange and no strange if they have fractal or integer dimension, respectively. 



In section 10 we will see how to compute the fractal dimension of an attractor. 
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(7„i37) solve (1/a = 



sqrt (2*%sigma/ (°/„beta* (%sigma +%rho) "2/2) ) , a) ; 



(%o37) [a = r+P±] 

2 / - 

v p 

(7„i38) solve(l/b = sqrt (2/ (%beta* (/.sigma +%rho) "2/2) ) ,b) ; 

(%o38) [b = t±A] 
2 - 

(7„i39) solve(l/c = sqrt (2*%beta/ (°/„beta* (%sigma +%rho) "2/2) ) , c) ; 

(%o39) [c = i^±4] 

These lengths depend only on the parameters <r, p and /3. Thus, taking R large enough 
(explicitly R> ^ max{2, \/T+7?, ^/l + ^/tr}), we obtain a sphere containing the initial 
conditions (here (x(0),y(0), z(0)) = (0,0,0) for simplicity) and such that when they evolve 
in time, all orbits remain inside this sphere. 



As already pointed, the presence of attractors having fractal structure (strange attrac- 
tors) is a signal for chaotic dynamics, but it is neither a necessary nor a sufficient condition. 
There are attractors derived from chaotic dynamics that are not strange (for example, in 
the logistic equation with k = 4 the attractor is the whole interval [0, 1]) and there are 
non-chaotic systems displaying attractors with fractal dimension (see [GOPY 84j ). 



7 Lyapunov exponents 

We have seen that a feature of chaotic dynamical systems is the sensitivity to initial con- 
ditions. Lyapunov exponents are introduced to qualitatively measure this property. The 
idea is quite simple: By choosing a fixed orbit, we compare it with other orbits having close 
initial conditions, and then measure the distance between them with a factor of the form 
exp(Ai). Here A is the Lyapunov exponent. If it is positive, the orbits diverge asymptot- 
ically and there is sensitivity to initial conditions (greater when when A is large). If the 
exponent A is negative, the orbits must asymptotically approximate each other and there 
is no sensitivity to initial conditions. A value of A = indicates that the orbit considered 
is a stable fixed point. 

Let us consider, for simplicity, a one- dimensional dynamical system 

Xn+l = f(x n ) 

with f : R — > M derivable except for a finite number of points (this is the case, for example, 
of the logistic map). We consider the orbits starting at the points xq and xq + e. After 
N iterations the points on the orbits will be, respectively, Jn(xo) and /at(xo + e), where 
In = f°- • •°f- So, after iterations the distance between the orbits is |/at(xo+£) — /at(xo)|, 
which can be written in the form 

\fN(xo + e)-f N (x )\=ee NX( ^ 
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for a suitable A(xo) 6 M. that is called the Lyapunov exponent at the point xq. Is in this way 
that the Lyapunov exponent gives a measure of how the initial separation e is amplified 
when the orbits evolve. 

Solving the last equation for A(xo): 

(%i40) assume (°/ epsilon>0,N>0)$ 
(°/ i41) remvalue(x[0])$ 

(7„i42) solve (abs (f [N] (x [0] +°/„epsilon) -f [N] (x [0] ) ) = 
%epsilon*exp(N*%lambda) ,%lambda) ; 



%o42) [A 



log 



|/jv(e+xo)-/iv(:Eo)| 



N 



and taking limits with e — > and N — > oo (taking into account the continuity of log and 
the definition of the derivative), we have: 



A(x ) 



lim — log 



df. 



N 



dx 



(*o) 



Just to show the symbolic capabilities of Maxima, we will use it to evaluate the derivative 
in the last expression. First, we declare the derivative of the function / by giving it a name 
(the usual /'): 

(5U43) gradef (f (x),f\ , (x))$ 

Then, we set the values x% = f(xo) and /'(xo) = ^(^o): 

(%i44) atvalue (f (y) , y=x [0] , x [1] ) $ 

(%i45) atvalue ( ' dif f (f (y) , y) , y=x [0] , f \ ' (x [0] ) ) $ 

Now, we compute the derivative of = f o f at xo (Maxima knows the chain rule): 
(%i46) at (dif f (f (f (x) ) , x) , x=x [0] ) ; 

(%o46) f (x ) /' (xi) 

It is then easy to prove (by induction) that 

df N 



N-l 



dx 



Oo) = Yl 



i=0 



so 



A(x ) 



lim — log 



N-l 

Yin 

i=0 



N-l 



lim — log I /' 

N^oo N ^ 1 



(7) 



i=0 



Using equation Q we can compute the Lyapunov exponents numerically, approximating 
the limit by a finite sum for iV big enough. As an example, consider the logistic map with 
parameter r = 3. In (%i8) we defined the function F(r, x), and now we define its derivative: 
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(7„i48) define (dF(r,x) ,diff (F(r,x) ,x))$ 



Then we set the parameter value: 
(°/.i49) r:3$ 

Here we set the number of iterations: 
(%i50) maxiter: 50000$ 

and the initial condition^ 
C/.i51) x[l] :random(1.0) ; 
(%o51) .4823905248516196 

Now, we construct the orbit from the initial condition: 
(%i52) for j:2 thru maxiter do (x[j] :F(r,x[j-l] ))$ 
Finally, we estimate the corresponding Lyapunov exponent: 
(°/ i53) L:0$ 

(%i54) for j:l thru maxiter do L:L+log(abs(dF(r,x[j] )))$ 
(%i55) '7olambda=L/maxiter ; 
(%o55) A = -3.145501884323275 1(T 4 

For this case, r = 3, we get A ~ 0. The reader can experimentally check that for the logistic 
map X(xq) = A, that is, the Lyapunov exponent does not depend on the initial condition 
xo G]0, 1[, it only depends on the parameter r. For values of r lesser than (approximately) 
3.569945, we have that A < 0, but if r > 3.569945, the behavior of A as function of r 
becomes quite complicated. 

8 The Duffing oscillator 

We have only considered systems described by first-order equations. However, most of the 
physical systems of interest are described by second order equations. In this section we 
present and example of these, the Duffing oscillator. 

Let us consider a general Newtonian system, described by a second-order differential 
equation: 

x(t) = F(x(t),x(t)) 

7 We denote it by xi instead of xq because we are going to make a list with the iterations, with the initial 
condition as the first element and, for Maxima, lists are enumerated starting with the index 1. 
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If F = F{x) = — dV/dx, we say that the system is conservative and the function V = V{x) 
is called the potential. For example, the system 

m = -\x\t) + X {t) 

is conservative and its potential is given by 
(%i56) 'V(x)=integrate(x~3/4-x,x) ; 



(%o56) V (x) 



A 1 
X X 



16 2 

The graph of this function is: 
(%i57) wxplot2d(integrate(x~3/4-x,x) , [x,-3,3] , [ylabel , "V(x) "] ) ; 




(%t57) 



%o57) 



We observe an unstable equilibrium at x = and two stable equilibria at x = — 2 and 
x = 2. Now let us modify the system by including a velocity-dependent damping term: 



x(t) 



1 



x 3 (t) + x(t) 



1 

To 



x(t) 



The system is no longer conservative, since the damping term dissipates energy in the form 
of heat. As a result, the oscillations made by the system decrease in amplitude until there is 
no motion. We illustrate this by computing the trajectory with initial conditions x(0) = 3, 
x(0) = 10: 

(%i58) duff : [v,-v/10+x-x"3/4]$ 
( /.i59) icduff : [3,10]$ 

(%i60) solduff :rk(duff , [x,v] , icduff , [t , 0, 100 , . 1] )$ 
(%i61) curveduf f : map (lambda ( [x] ,rest(x,-l)) , solduff )$ 
(%i62) pointsduff : map (lambda ( [x] ,rest(x)) , solduff )$ 
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(%i63) wxdraw2d(point_type=none , point s_ jo ined=true , 
xlabel="t" ,ylabel="x(t) " , points (curveduff) ) ; 




20 40 60 80 100 



(%t63) * (%o63) 

We can also represent the dynamics of the system in the plane (x, v = dx/dt) (the phase 
plane): 

(%i64) wxdraw2d(point_type=none,points_joined=true , 

xlabel="x(t) " ,ylabel="v(t) " .points (pointsduff ) ) ; 




Note that this result coincides with (%o63), where it is clear that the system evolves to the 
stable equilibrium x = —2, oscillating around this value with decreasing amplitude until it 
stops. 

The Duffing oscillator is obtained including in the previous system and oscillating force 
(the forcing term) sin(wi): 

1 , 1 

x(t) = --x d (i) + x(t) - j^x{t) + sin(cji) 

The effect of the forcing term is to restore the energy dissipated by the damping term. 
If the frequency of the forcing term is the appropriated one, the system will oscillate in 
a stable way, but if there is no synchronization with respect to the non-forced system, a 
chaotic behavior will appear. For example, for oj = 1 starting from the unstable equilibrium 
(x = 0, v = 0), after a short transient the system will evolve to a stable oscillating regime 
(a limit cycle): 
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(7,i65) duffing: [v,-v/10+x-x~3/4+sin(t)] $ 

(7„i66) icduffing: [0,0]$ 

C/.i67) sduffing:rk (duffing, [x,v] , icduffing, [t , 0, 100 , . 1] )$ 

(%i68) cduff ing: map (lambda ( [x] ,rest(x,-l)) ,sduff ing)$ 

(%i69) pduf f ing: map (lambda ( [x] ,rest(x)) ,sduffing)$ 

(%i70) wxdraw2d(point_type=none , point s_ jo ined=true , color=magenta, 
xlabel="t" ,ylabel="x(t) " , points (cduff ing) ) ; 



(%t70) 




20 40 60 

t 



80 100 



%o70) 



The corresponding phase plane diagram shows that the system tends to a limit cycle, 
oscillating around the two stable equilibria, instead of approaching just one of them, as 
happens for the non- forced case: 

(%i71) wxdraw2d(point_type=none,points_joined=true , color=magenta, 
xlabel="x(t) " ,ylabel="v(t) " .points (pduf f ing) ) ; 




The behaviors considered do not exhaust the possibilities for the dynamics of the Duffing 
oscillator. In the next section we will see how to analyze the chaotic case. 
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9 Poincare sections 



In general, for a second order differential equation of the form 

x(t) = F(t, x, x) 

or, equivalently (introducing the variables u(t) = x(t), v(t) = x(t)) for the first order system 

v 

F{t,u,v) 

the trajectories in the phase plane can intersect themselves and the resulting diagram be- 
comes quite complicated. Let us note that this is not a contradiction with the theorem on 
the uniqueness of solutions, because here we have a non autonomous differential equation 
and the phase plane is obtained projecting from the phase space (t,x,v) h-» (x, v). The 
theorem applies, when it is the case, in the variables (t,x,v). For example, for the Duffing 
oscillator with forcing term 2.5sin(2t), we get the following phase diagram: 

(%i72) duffingl: [v,-v/10+x-x~3/4+2.5*sin(2*t)] $ 
C/.i73) icduffingl: [0,0]$ 

(°/ i74) sduffinglrrk (duffingl, [x,v] , icduffingl, [t ,0 , 100,0 . 1] )$ 

(%i75) pduf f ingl :map(lambda( [x] ,rest(x)) ,sduff ingl)$ 

(%i76) wxdraw2d(point_type=none , point s_j oined=true , color=coral , 
xlabel="x(t) " ,ylabel="v(t) " .points (pduf f ingl) ) ; 



' du 
dt 

dv 
- ch 7 




Poincare introduced a very useful technique to study the dynamics in these situations. 
It consists in reckoning a trajectory in the phase space (x, v, t) by using hyperplanes equally 
spaced in time as "check points", and then projecting the resulting points on the plane 

(x,v): 

Of course, if the time interval between these hyperplanes (called the Poincare strobo- 
scopic sections) is T and the system is periodic with period T, we will see just a point on the 
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Figure 1: Poincare sections 



phase plane. In the figure we see a quasi-periodic system: on the phase plane the trajec- 
tories are circumferences, but several points can appear in the sections projection, showing 
that the system oscillates between several equilibria with a frequency different from the 
time interval between the Poincare sections. 

When the dynamics is chaotic, the diagram obtained by projecting the Poincare sections 
can be really complicated, usually displaying fractal structure. 

Let us compute the Poincare sections for the case of the Duffing oscillator. First, we consider 
the case of a period T = ~ corresponding to the frequency u) = 1 (that is, the forcing term 
is sin(t)). We register the points in the orbit at time intervals of length T = — = 2tt: 

C/.i77) miter: 25$ 

(°/.i78) 7„tau:bfloat(7„pi)$ 

(70179) sduf fing2 : rk (duff ing, [x,v] , icduffing, 
[t , ,miter*2*7otau , 7.tau/30] ) $ 

(70180) pduf f ing2 : create_list (sduf f ing2 [i] , i ,makelist (i*60 , i , 1 ,miter) ) $ 

(70181) poinduf fing2 : map (lambda ( [x] ,rest(x)), 
makelist (pduf f ing2 [i] , i , 1 , miter) ) $ 

(70182) wxdraw2d(point_type=f illed_circle , color=magenta,xtics=l , 
ytics=l ,xrange= [-4, 1 . 5] ,yrange= [0,4] .points (poinduf fing2) ) ; 
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I , , , , , 1 

4-39101 

(%t82) (%o82) 



With the exception of some sparse points representing the transient, we see that all the 
points tend to a cluster (oscillations occur close to a limit cycle). To show that the most 
distant points corresponds to the transient, we eliminate, for example, the first four samples 
(corresponding to the sections at t = 0, 2tt, 4tt and 6tt): 



(%i83) transient : map (lambda ( [x] ,rest(x)), 
makelist (pduf f ing2 [i] , i , 5 , miter) ) $ 



(%i84) wxdraw2d(point_type=f illed_circle , color=magenta,xtics=l , 
ytics=l ,xrange= [-4, 1 . 5] ,yrange= [0,4] .points (transient) ) ; 
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1 



I 1 1 1 . . 1 

4-39101 

(%t84) (%o84) 

Of course, increasing the number of samples we get a more illustrative picture. But the 
computational effort is not worth in this case. There are more interesting situations, as it is 
the case of the Duffing oscillator with forcing term 2.5sin(2i). Now the period corresponds 
to the frequency uj = 2, that i S , T = \ = 7T. We let the system evolve for a long time in 
order to get a large sample via the Poincare sections: 

(°/ i85) duffing3: [v,-v/10+x-x~3/4+2.5*sin(2*t)] $ 
C/.i86) icduffing3: [0,0]$ 
(%i87) maxiter : 1000$ 
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(%i88) sduff ing3:rk(duff ing3, [x,v] , icduf f ing3 , 
[t , ,maxiter*°/ tau , %tau/30] ) $ 

( /.i89) pduffing3: 

create_list (sduf f ing3 [i] , i ,makelist (i*30 , i , 1 ,maxiter) ) $ 

(%i90) poinduf fing3 : map (lambda ( [x] ,rest(x)) , 
makelist (pduf f ing3 [i] , i , 1 ,maxiter) )$ 

(%i91) wxdraw2d(point_size=0 . 3 ,point_type=circle , 

xrange= [-5 , 5] , yrange= [-7 , 3] , xtics=l , ytics=l , 
color=coral .points (poinduf fing3) ,grid=true) ; 



(%t91) 



% **t V* ■ 

■ ■ -a* Jf* : t*> 



-: # - -:- *+f& 
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k:. 

: : • : : . i? •*: : 



oo91) 



10 Fractal dimension 

As we have pointed out in previous sections, the appearance of strange attractors is a usual 
feature of chaotic systems. These are sets around which the orbits move asymptotically 
(thus justifying the name of attractor) and frequently have fractal dimension. We will 
introduce a particular notion of fractal dimension (there are other possibilities) in the case 
of a two dimensional dynamical system, so the attracting set A will be a subset of the plane. 
An example of this situation is given by the Duffing oscillator. The box-counting dimension 
of A is defined as follows: we start with a square of side length / enclosing A. In the next 
step, we divide the square in r sub-squares, each having side length l/r, and so on. In the 
&-step we have a grid with r k ~ l squares, each with side length 5k = -jc=i: that cover the 
set A as in (%o91). For each step k, let N(k) be the number of squares containing at least 
some point of A. 

Let us observe that if A were a smooth curve (and, therefore, homeomorphic to a segment 
of R, to which we would assign a topological dimension 1), for a small enough value of the 
side of the squares each of them would contain a piece of the curve with length as close as 
we wish to the side length of the square. So, in the case of a smooth curve, we would have 

lim N(k) -S k = L, 

k— >oo 



22 



with L the length of the curve. In this case we say that the scale of N(k) goes as 67 , 
that is, N(k) ~ 87 . Analogously, if A were a measurable set on the plane (a set having 
an area and, thus, of dimension 2), we would have N(k) ~ 57 2 . Let us note that in these 
cases, that can be called "regular" ones, the exponent D in the relation N(k) ~ 57 D equals 
the topological dimension of the set. Then, a fractal set can be defined as a set for which 
a relation of the kind N(k) ~ 57 D holds, with D a positive real number (not necessarily 
an integer). In general, we call this number D the box counting dimension of the set, and 
remar that it satisfies 

k^oo — log Ofc 



This is a constructive definition that suggest how to perform the calculation of the box- 
counting dimension for a given set. Explicitly, we can plot logiV(A;) versus — log^fc, then we 
can fit the data and finally estimate D as the slope of the regression line. As an example, 
let us calculate the dimension of the Duffing attractor (%i85) with this procedure. We note 
that the Duffing attractor is contained in the box [—5,5] x [—7,3]. We divide this box in 
10 x 10, 20 x 20,..., 200 x 200 boxes, counting how many points are inside them. As a detail, 
we first "normalize" the coordinates in such a way that they lie in the box [0, 10] x [0, 10]. 
The algorithm starts with a matrix filled with zeros and changes the scale at successive 
steps multiplying these lengths by 2, 3, 20: 

(%i92) resolution: 20$ 

(%i93) for n:l thru resolution do 

(Z[n] :substpart(" [" ,zeromatrix(10*n, 10*n) ,0) , 
boxcount [n] : , 
for k:l thru maxiter do ( 
ix : floor (n*(poinduff ing3 [k] [l]+5)) , 
iy : floor (n*(poinduff ing3 [k] [2] +7)) , 
if is(Z[n] [ix] [iy]=0) then 
(Z[n] [ix] [iy] :1, boxcount [n] : boxcount [n] +1) ) )$ 

Here is the plot of the pairs of numbers obtained. 

(%i94) f itdim:makelist( 

[log (n) /log (10) , log (boxcount [n] )] ,n, 1 .resolution) ,numer$ 

(%i95) wxdraw2d(point_size=l ,point_type=f illed_circle , 
color=dark_violet , points (f itdim) ) ; 
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(%t95) 




(%o95) 



The first values are "outliers" (since the boxes are in fact too big), so we can safely discard 
them. In the same way, for the last values the boxes are too small, counting nearly one 
point per box, so again we can discard them. This situation, of data that gives a wrong 
contribution to the estimation of the box-counting dimension, is a well known problem 
discussed in the literature (see |Kli 9 4]); in fact, there are several algorithms aimed to 
automatically selecting the most representative data, but for our (limited) purposes, we 
just use the visual information in the graph (%o95). According with this, it seems a good 
choice to eliminate the first 7 values and the two last ones to make the fitting. The result 
is: 

(%i96) f itdimension:rest(rest(fitdim,7) ,-2)$ 
(%i97) load(stats)$ 

(%i98) simple_linear_regression(f itdimension) $ 

(%i99) f iteqn:take_inf erence( 

model , simple_linear_regression(f itdimension) ) ; 

(%o99) .6534032420835236 x + 6.021481801936183 
C/oilOO) fractaldimension: coef f (f iteqn.x) ; 

(%ol00) .6534032420835236 



Since our approximation is very rough, we can round-off the value so obtained. Therefore, 
we can say that the fractal dimension of the Duffing attractor is approximately 0.65 (between 
and 1, as it should be expected). 
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